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Using a collective Monte Carlo algorithm we study the low-temperature and long-distance prop- 
erties of two systems of two-dimensional classical tops. Both systems have the same spin-wave 
dynamics (low-temperature behavior) as a large class of Heisenberg frustrated spin systems. They 
are constructed so that to differ only by their topological properties. The spin-stiffnesses for the 
two systems of tops are calculated for different temperatures and different sizes of the sample. This 
allows to investigate the role of topological defects in frustrated spin systems. Comparisons with 
■ Renormalization Group results based on a Non Linear Sigma model approach and with the predic- 

tions of some simple phenomenological model taking into account the topological excitations are 
done. 
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to ; T. INTRODUCTION 

+-* 

S The long-distance behavior of the classical Heisenberg AntiFerromagnet on a Triangular lattice (HAFT model) has 
been the subject of nauch interest. In three dimensions a most important issue is the nature of the universality of 
*^ | its phase transition uu . In two dimensions, this model has also been widely studied since it exhibits a non-trivial 
finite-temperature behavior due to the presence of topological excitations. Topology enters the problem since the 
order parameter of the model belongs to SO(3) whose first homotopy group is ttx[SO(3)] = Z^. As a consequence, 
there exist topologically stable point defects - called vortices - for this two-dimensional system. Arguments involving 
entropy and energy of the defects suggest the occurence of a change of behavior at a finite temperature Ty between 
£N) , a pure spin- wave regime with confinea^-vortices for T < Ty and a regime of free vortices for T > Ty. Several Moate 
£> ■ Carlo studies of the HAFT model W~t3 and of some generalizations with an easy-axis exchange anisotropy GdiU 
C^ 1 | have indeed revealed the existence of various regimes resulting from the presence of defects. 

Here, our purpose is to shed some light on the interplay between vortices and spin-waves in 2D by studying 
with Monte Carlo simulations two lattice models of ferromagnetically interacting tops. Both models have the same 
spin-wave dynamics as the original HAFT model but they differ by their topological properties: the first one has 
the same topological content as the HAFT model, the second one is topologically trivial. The role played by the 
topological defects emerges from the comparison between these two models. Note that this comparative study would 
have been more difficult to implement directly on the original HAFT model. The physical quantity we consider 
in our study is the spin-stiffness which, for a spin system jMft-f finite lattice of size L, measures the free-energy 
increment resulting from a twist of the boundary conditions EaO . The spin- wave part of the spin-stiffness identifies 
with the coupling constant of the Non Linear Sigma (NLcr) model renormalized at scale L by thermal fluctuations. 
Accordingly, the behavior of the spin-stiffness as a function of the size of the lattice provides a direct test of the 
perturbative Renormalization Group (RG) predictions of the NLcr model. For the models of tops studied here we 
expect that, at sufficiently low temperature -for T significantly smaller than Ty- where the physics is dominated 
by pure spin-waves, the behavior of the spin-stiffnesses agrees with the RG predictions. On the other hand, we also 
expect that, for a topologically non-trivial model, near Ty, the vortices also contribute to the spin-stiffness. For such 
a model the behavior of the spin-stiffness must disagree with the standard RG predictions. One great advantage 
of the spin-stiffness p is that, in contrast with the correlation length £ which cannot be easily computed at very 
low temperature since it diverges typically as exp(l/T), p has a smooth behavior at low temperature. It is thus, in 
principle, easily computable. Regarding Monte Carlo simulations, a central aspect is that, at the low temperatures 
we are interested in, the dynamics of 2D spin systems is governed by strongly correlated spin-waves, independently 
of the presence of vortices. These modes are responsible for a severe critical slowing down which makes difficult the 
convergence of simulations based on local algorithms in which one spin is flipped at each Monte Carlo step ("local 
update" Monte Carlo-schemes). To resort to collective algorithms based on global updates (construction of clusters) 
is then important EDO. However, as well-known, such algorithms work well for ferromagnetic systems but not 
for frustrated ones. Note that, together with the ability of comparing topologically different models, this aspect is 
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an additional motivation to consider ferromagnetic top models rather than the original HAFT one. Actually, the 
implementation of the basic rules of collective algorithms for systems consisting of tops is itself not so clear. However, 
this difficulty can be circumvented by rewriting the models of tops considered here as ferromagnetic four-component 
spin systems while, of course, preserving both their spin-wave and topological contents. Thanks to these various 
tricks and to the cluster algorithm we are then able to scan a large temperature range below Ty while fully controlling 
the convergence of our simulations. It should be noted that previous Monte Carlo calculations of the spin-stiffness 
on related .models (HAFT model and generalizations) have been done by using Monte Carlo schemes with local 
updates iT~Ej. Note that, in contrast with these works, we have considered here very small temperatures, T <C Ty. 
To resort to a non-local algorithm to accelerate the convergence of simulations is therefore essential. Note also that 
at the intermediate temperatures where the necessity of using global Monte Carlo schemes is less important we have 
also found a clear improvement associated with the use of cluster algorithms. 

The main result of this paper is that we have found some striking differences in the behavior of the spin-stiffnesses 
as a function of the linear size for the two models, with and without topological excitations. At very low temperatures 
the temperature-rescaled spin-stiffness, p = p/T (the natural quantity to consider, see below) of both models displays 
the characteristic behavior: 

predicted by the perturbative RG approach of the NLcr model @. This is, of course, expected since the spin- wave 
contents of both models are identical. We call the regime corresponding to this range of temperatures the "spin-wave" 
regime. Note that, for one value of the temperature, this asymptotic scaling of p with respect to InL has already 
been confirmed directly on the HAFT model by Southern and Young . At higher, but still low, temperatures the 
two models begin to display different behaviors. While the spin-stiffness of the topologically non-trivial model still 
displays the previous characteristic behavior, its absolute magnitude with respect to the trivial model is found to 
decrease quite rapidly as a function of the temperature. We propose to refer to this regime as an "almost-spin-wave" 
regime, a regime where the only significant effect of vortices is just to shift down the value of p/T. Next, at higher 
temperatures we enter a regime called here the "vortex" regime where the vortices play a major role. In this regime, 
the spin-stiffness looses its regular behavior. It exhibits large fluctuations around its mean value with the presence 
of "plateaux" and abrupt jumps as a function of the linear size. Nevertheless, by considering the global behavior of 
the curve it is still possible to define some effective linear regime as a function of In L similar to that described by 
Eq. ([!]). However, in contrast with the spin- wave regime, the slope of the spin-stiffness is no longer constant (l/47r) 
and is found to increase quite rapidly as a function of the temperature. It is remarkable that this regime is observable 
only within a narrow range of temperatures. At slightly higher temperatures, the curve of the spin-stiffness recovers a 
much more conventional behavior: smooth decrease as a function of the size and cancellation of p at some finite lattice 
size corresponding to some finite correlation length. Regarding the theoretical interpretation of our results, we show 
that the very low-temperature regime is in full agreement with the RG predictions. The so-called "almost-spin-wave" 
and "vortex" regimes are much more puzzling. However, it is shown that the most salient features induced by the 
topological defects in these regimes can be rather well reproduced using some simple phenomenological RG equations 
which combine the topology of the XY model and the spin-wave content of the 0(4) model. 

The organization of the paper is as follows. In Section II, the various actions of the lattice models are presented. 
In Section III, the expressions of the spin-stiffnesses suitable for Monte Carlo simulations are given. In Section IV, 
we present briefly the Wolff-Swendsen-Wang algorithm used. Our results are given in Section V. In this latter section 
we present the behavior of the spin-stiffnesses as a function of the lattice size in the various temperature regimes 
going from low- to high-temperatures. Finally, we present in the last section our first attempt toward a theoretical 
interpretation of the effect of the vortices in the almost-spin-wave and vortex regimes. 



II. THE LATTICE MODELS 

A. The SO(3) <8> 0(2) top model 

Our first step is to,map the HAFT model into an equivalent non-frustrated one. As shown by Dombre and Read 
and Azaria et al. O the long-distance effective hamiltonian of the HAFT model consists in a system of classical 
interacting tops. This can be understood from the fact that the 120° structure of the spins of the HAFT model in 
the ground state fully breaks the SO(i) symmetry so that the order parameter is a rotation matrix R £ SO(3), a 
classical top. As in the non-frustrated case, once the theory is reformulated in terms of ±]ie_order parameter, the 
effective interaction becomes ferromagnetic. The hamiltonian of the top model thus reads EjEJ : 
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Hi = - Yl Tr (PKi'Rj) (2) 

<i,j> 

where Ri is a rotation matrix of 50(3) defined on site i and P=diag(p 1 ,p 1 ,p 3 ) is a diagonal matrix of positive 
coupling constants which represents the interaction strengths between the different axes of the topfcpNote that the 
temperature has been included in the piS. The HAFT model corresponds to the special case»3 = . The 50(3) 
symmetry of the HAFT is realized here through the rotational invariance of hamiltonian (||) under left global 5(3(3) 
rotations Ri — > URi, U G 50(3). With the matrix P considered here it is also invariant under the 0(2) group of 
right global transformations: Ri — > RiV that commute with the matrix P. Thus, hamiltonian (|J) is invariant under 
the group G = 50(3) (g> 0(2). This left 0(2) group is reminiscent of the C% v symmetry of the triangular lattice. Note 
that it will be convenient in the following to consider the case p 3 ^ sincapthe hamiltonian made with this P is the 
general one invariant under G = 50(3) ® 0(2) and that, as well known the condition p 3 = is not preserved 

by renormalization. 

The symmetry breaking pattern described by hamiltonian (Q) is G = 50(3) (g> 0(2) broken down to H = 0(2): 

o = so^m sSm (3) 

where the notation = means that G/H is topologically isomorphic to 50(3). The symmetry breaking pattern 
thus corresponds to a fully broken 50(3) group. For the original HAFT model, this symmetry breaking pattern is 
G/H = 50(3) (g) C^vlCzv and is thus identical to that given by Eq.(|J). This is the reason why the substitution of 
the discrete C^ v by the continuous 0(2) one in hamiltonian (|2|) is harmless. It is interesting to note the identity: 
50(3) = 50(4)/(50(3) ® Z2) = S3/Z2, S3 being the 3-sphere, since it shows that the model of tops (Q) is equivalent 
to that of four-component spins living on the four-dimensional unit sphere with antipodal points identified. This will 
allow us, in the following, to build a vector model equivalent to the preceding matrix one and suitable for Monte 
Carlo simulations. Note finally that when p\ = p 3 the symmetry group is enlarged to G — 50(3) <E> 50(3) and the 
symmetry breaking pattern is G/H = 50(3) <8> 50(3)/50(3). 

B. The 5(7(2) ® U(l) top model 

We now build the topologically trivial counterpart of the previous top model. We want to preserve the spin-wave 
part of the model while discarding the topological excitations. Since the spin-wave excitations correspond to small 
fluctuations of the order parameter, they only probe the local structure of the order parameter space G/H and not 
its globaL|p^topological - structure. This local structure is itself completely determined by the Lie algebras of G 
and H HnJ. We thus need a model defined by an order parameter space G 1 / H' locally isomorphic to G/H and 
topologically trivial. This is obtained by considering the covering group 5(7(2) of 50(3). The relevant model is thus 
built on the manifold 5(7(2) (g) (7(1)/(7(1). The most general hamiltonian invariant under SU(2) ® (7(1) writes: 

#2 = - J2 l 2 ^ 1 + Tr 9^9j + \(Pi- Ps) (Tr a 3 g^g.fX (4) 
<i,j> J 

where <?.; € SU (2) and 03 is the third Pauli matrix. The first term in this hamiltonian is clearly invariant under the 
cross product of a left 5(7(2) group and a right 5(7(2) group: gi — > MgiN, M, N S 5(7(2). The second one explicitly 
breaks the right 5(7 (2) down to a right U (1) so that the hamiltonian is generically 5(7(2) ® U (1) invariant. This (7(1) 
symmetry corresponds to the 0(2) symmetry of hamiltonian (E|). 

The symmetry breaking pattern described by hamiltonian (Q) is: G' — 5(7(2) ® (7(1) broken down to H' — (7(1): 

so that it corresponds to a fully broken 5(7(2) group. Again, it is interesting for the following to note the identity 
5(7(2) = 50(4)/50(3) = S3 which means that the model (|j) is equivalent to that of four-component spins living on 
the four-dimensional unit sphere. Again, when p\ = ^3, the symmetry group is enlarged to G' — 5(7(2) ® 5(7(2) and 
the symmetry breaking pattern becomes 5(7(2) £g> 5(7(2)/5(7(2). Note finally that the choice of coupling constants in 
(|J) is such that the two models (0) and (0) have the same temperature scale. 

The models corresponding to (0) and (B) have, by construction, the same spin-wave dynamics but can neverthe- 
less strongly differ when excitations associated with the topology are activated. The 5(7(2) ® (7(1) model being 
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topologically trivial, i.e. iri(SU(2)) = 0, we expect that a Monte Carlo study of this model will be well reproduced 
by a pure spin-wave approach. We show in the following that this is indeed what happens: as in the topologically 
trivial 0(N)/0(N — 1) ferromagnetic spin systems, the critical properties of the SU(2) Cgi i7(l) model are in perfect 
agreement with the perturbative RG predictions made on the continuous limit of the top model, a NLu model. On 
the other hand, the 50(3) (g) 0(2) model being topologically non-trivial, some disagreement between the perturbative 
and Monte Carlo approaches at sufficiently high temperatures are found, as expected. 



C. The vectorial version of the SO(3) <g> 0(2) model 



As already mentioned in the introduction, the cluster algorithms are easier to implement for spins than for matrices. 
We thus need vectorial versions of our hamiltonians. This is achieved by using the decomposition of a rotation matrix 
Ri of 50(3) in terms of a four-component unit vector 5 t = (£?, S l ) = (Sf, S},Sf,Sf): 



The hamiltonian 



rykl n / ok qI 

then takes the form 



<hj> 



n-i 



2eki m SfS™ +2(5'. 



(6) 



+ 4(p 3 - pi) 



(s®Sj + sfSj) (slsj + sfSj) + (sfSj — SjSf) (sis] — sfsj) 



(7) 



QO I c3 2 cl 



si 



sf 



s) 



The first term of (^) represents the hamiltonian of a system of spherical tops - i.e. SO (3)® SO(3) ~ SO(4) symmetric 
- for which p± = p%. This term is also invariant under a local - gauge - Z2 group. This Z2 symmetry expresses the 
non trivial topological character of the 5*0(3) group. The first term of ([jj) is also known as the hamiltonian of the 
RP 3 — SO(4)/ (SO(3) <g> Z2) model which expresses the isomorphism between the manifolds 50(3) and RP 3 . This 
model and, more generally, the jS/pn models for general N have been extensively studied and the question of the 
nature of their continuum limit IE3~E3 strongly debated, also in connection with topological defects. The second term 
of hamiltonian (0) also displays the Z2 local symmetry but breaks the global 50(4) symmetry so that the hamiltonian 
is generically globally 50(3) ® 0(2) and locally Z 2 symmetric. 



D. The vectorial version of the SU(2) ® U(l) spin model 

It is also possible to express hamiltonian (Q) in terms of four-component vectors Si- Using the decomposition of a 
SU(2) matrix: 

9i = S° + ia.Si (8) 
ct/c, k — 1, 2, 3 being the Pauli matrices, hamiltonian (|4|) writes: 

H 2 = - E 4 (Pi + J*) + 2(pi - p 3 ) (S°S 3 ~ S°S 3 - SfSj + S]Sl) 2 . (9) 

<i,j> 

In this expression the first term is 0(4) globally invariant and corresponds to the hamiltonian of a four-component 
ferromagnet whereas the second term breaks this symmetry. The hamiltonian (^|) is thus generically globally SU (2) g) 
U(l) symmetric. Note that the Z2 local symmetry has now disappeared. This is a consequence of the trivial 
topological character of SU(2). Note also that the scale of temperature has been chosen so that both Hamiltonians 
H[, Eq.(|), and H' 2 , Eq.(@), have the same linearized spin- wave form in the symmetric case (pi = p$), namely 
PI' X I2 — ~ ^2<ij> SpidSi.SSj where SS represents the spin deviation from the reference vector. Finally, remark that 
the Jacobians resulting from the change of variables: matrices — > spins in both 50(3) (g> 0(2) and SU(2) U(l) are 
trivial and thus do not contribute to the free energy. 
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III. THE SPIN-STIFFNESSES 



A. The spin-stiffnesses of the lattice models 

The spin-stiffness p a to be computed numerically is defined as the free energy increment under twisting the boundary 
conditions, for instance in the x direction around the direction a. This is realized by coupling the system with two 
walls of tops: R(x = 0) = R\ and R(x = L) = R2, R2 being deduced from i?i by a rotation of angle 9 a around the 
direction a and by measuring the variation of the free energy with respect to 9 a : 



d 2 F(9 a ) 



For a system with partition function: 



we have: 



-T 



[Ri] 



d 2 H 



H 



I dH\ 2 

~\doJ \oei 



f dH 

\de~ a 



(10) 



(11) 



(12) 



o a =o 



Since H is even in the (9 Q 's, the average value of dH/d6 a is equal to zero and only the two last terms of Eq.jT^) 
need to be computed. 

1. The spin-stiffnesses 0} the SO(S) ® 0(2) model. 
In principle we have to compute the different average values in (n3) from the partition function: 



[Ri] <t,j> 



constrained by the boundary conditions: 



R(x = 0) = i? 

R{x = L) = R a e i9aTa 



(13) 



(14) 



where Rq is a rotation matrix of reference (e.g. Rq = I), T a is the generator of rotation around the a direction, and 
9 a the angle of rotation (Eq.([l4|) must be understood without the sum over a). However, in practice, the presence of 
derivatives with respect to 9 a in Eq.(|l2|) as well as the fact that the cluster algorithm is implemented with spins, makes 
expression ( |l2| ) not suitable for our simulations. We proceed in two steps to reformulate the model in a numerically 
convenient way. First, to get rid of the derivatives, we compute them analytically and rewrite the average values in 
( |l2| ) as ^-independent quantities. To do this we decompose Ri into a zero temperature part - Rf - and a fluctuation 
part - hi - : 



Ri — Ri hi 



(15) 



where both Rf and hi belong to 5*0(3). In Eq. (|l5|) , Rf is by definition a solution of the classical equations of motion 
and thus reads: 



Rf=e L 

and the hi's satisfy the boundary conditions: 

h(x = 0) = h(x = L) = 1. 

We thus have: 



(16) 
(17) 



5 



dH 



= \ E Pk e almhi l h^ m (Xi - Xj) . 
8a=0 <i,j>,k,l,m 



(18) 

~ <i,i> ,k,l,m 



The average values in Eq. ( |l2] ) must now be computed with 



Z = E ex P( E T^V 1 ^) (19) 

[hi] ^ <i,j> ' 

with the boundary conditions (|l7|). At this stage, everything is written in terms of the h's. Thus we can now perform 
the second step of our derivation that consists in using in Eq.®, (jlf) and (|l|) the same decomposition as in Eq.ra) 
but now for the h^s: 

hf = 2 (s$Sl - \s kl ) + 2e kln S°Sr + 2 (V - \) hi • (20) 



Since now all thermal average values are entirely expressed in terms of four-component spins, the cluster algorithm 
can be implemented to compute the spin-stiffnesses. 

2. The spin-stiffnesses of the SU(2) (g> U(l) model. 
The same method can be employed for the SU(2) (g) U(l) model for which we have: 

Z = E cx p( E 2 (P! + Tr g-'gj + i(p x - p 3 ) (Tr g^gtf ) (21) 

[Si] V <i J> J 

with, again, the fixed boundary conditions: 

g(x = 0) = g 

%bJ± (22) 

5 (x = L) = g e 2 

where go is a rotation matrix of reference of SU(2) (e.g. go = 1), is a Pauli matrix, and 8 a the angle of rotation 
(Again in Eq.(^) there is no sum over a). 

As in the SO(3) case, we make the decomposition in classical and fluctuating parts: 

9i = gf K (23) 

with: 

. n &a Xi 

gf = e 2 L (24) 

and h satisfying: 

h(x = 0) = h(x = L) = 1 , (25) 



The different terms of (|12|) are separated into 50(4) and SU(2) ® £7(1) symmetric parts. Writing S — Sso(4) + 
Ssu(2)®U(i) an d using the decomposition: 

fH = S° + ia.Si , (26) 

we have with obvious notations: 
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(P3 



0(4) 



s a =o 



oE 



PlJ ^ ~y 

<ij> 



Sj.Sj 



P3 



<ij> 



Si-Sj 



S4 S j 



S}S S 



C2 CO CO Cl 
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) E 34- 



(Xj Xi) 



(27) 



with: 



SU{2)®U(1) 



*=0 



P3J 



■ E 

<*,j> 



ol c3 c2c0 c0q2 



\ (SiS® + SfSj — S i Sj — S i S, S[ 



qO q3 



c3 qO ni qz 1 q'z q 
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The spin-stiffnesses can now be computed with: 



E exp ( E 2(pi +p 3 ) Tr ftr 1 ^ + i(p x -p 3 ) (Tr <7 3 frr 1 /^) 2 ) 



with the set of equations ( |27| ) and the boundary conditions fl25|) or, in terms of the four-component vector S: 

S(x = 0) = 5(a; = L) = (1, 0, 0, 0) . 



(28) 



(29) 



(30) 



B. The spin-stiffnesses of the NLct models 



The spin-stiffnesses can be analytically computed from the continuum versions of the top models, which is a 
Non Linear Sigma (NLcr) model. They identify with the effective coupling constants of the NLcr model at scale L, 
renormalizcd by thermal fluctuations. Let us recall that the perturbative treatment of the NLcr model takes only 
into account the spin-wave part of the spin-stiffnesses. Indeed, the (3 functions of any NLcr model are completely 
determined by the local properties, i.e. by the metric, of the manifold G/H while they are insensitive to its global, 
i.e. topological, structure. Since the metric itself is completely determined by the Lie Algebras of G and H, the 
perturbative (3 functions - and thus the behavior of the spin-wave part of the spin-stiffnesses - of the| r Sp(3) <£> 0(2) 
and SU(2) Cg) U(l) NLcr models are identical by construction. They are given at two-loop order by: 1 



dpiil) 
dl 

dm) 



1 



1 Ml) 5 p 3 (l) , 3 p 3 il) 



1 



2tt 4tt?i(Z) 327r 2 p 1 (?) 3 8^ 2 ^{if ^ 2 Pi{l) 



1 p 3 (lf 1 p 3 (lf 



(31) 



dl ^p!(l) 2 32 ^Pl(0 4 

with I = hiL/a where L is the system size, a the lattice spacing, the initial conditions of the RG flow being given by: 



pi(l = 0) = pi +p 3 
p 3 {l - 0) = 2 Pl . 



(32) 



Note that the RG equations ( pl| ) arc written in terms of quantities that contain, in their definition, the temperature, 
as it is clear from Eq. ( |32"| ) since the temperature is included in the pi 's. The Monte Carlo counterpart of the p a 's are 
thus given by the p Q /T's previously defined. 



7 



In contrast with the O(N) case, there are a priori three different spin-stiffnesses in our case but, in fact, only two 
are independent since the left 0(2) symmetry constrains two of them to be identical. Note that, asymptotically, i.e. 
for L ^> a, the behavior of the p a 's is given by the infrared limit of the flow equations. It is easy to show that, in this 
limit pi — + P3 and the WfxJel becomes effectively SO(3) ® SO(3) ~ SO (A) symmetric. This is the known phenomenon 
of enlarged symmetry Q>E!l . We recover therefore the universal scaling of the spin-stiffness of a ferromagnetic N = 4 
vector model. At leading order we have: E3 

1 £ 

pi ^p 3 ^_ m _ (33) 



where £ is the correlation length. Note, of course, that since the spin-stiffnesses in Eq.(31 ) are calculated perturbatively, 
they can only be valid in a limited range of low temperatures where the perturbation theory at two-loop order is 
meaningful. 

From the preceding analysis, it should be clear that for a topologically trivial model, the spin-wave part of the spin- 
stiffness p a identifies with p a /T, so that we can expect that this last quantity follows the RG equations calculated by 
perturbation theory. Indeed, for the O(N) model, it has been checked that, at sufficiently low temperature^, the 
spin-stiffness p a /T follows the RG equations calculated by means of the 0(N) NLc model up to two loop order ESta. 
In the same way, we expect the different spin-stiffnesses p a /T of the lattice SU(2) ® U(l) model calculated by Monte 
Carlo simulation to follow the RG equations (|3l|) in a large range of low temperatures since it is topologically trivial. 
On the other hand, we expect the behavior of the spin-stiffnesses p a /T of the lattice SO(3) 0(2) model to agree 
with ( |3"l| ) at very low temperature, i.e. below TV, where the topological defects are not activated, but to disagree 
with the perturbative RG predictions near and above the cross-over temperature Ty. 



IV. WOLFF-SWENDSEN-WANG ALGORITHM 

The simulations presented in this paper are based on a generalization of the Wolff-Swendsen-Wang JlTj , |l8| . [Q , 
pl| algorithm to N- vector models as presented by Caracciolo et al. in j32|. The method is based on an embedding 
of Ising spins e into the iV-component (here N = 4) continuous spins S according to: 

S t = S^ + e,\sf\r (34) 

where r is a unit vector chosen randomly on the sphere S 3 , = Si — (S{ • r)f and Sf — (Si -r)f are the components of 
the spin vector perpendicular and parallel to the unit vector r, the Ising variable being given by ti — sgn(Si-r) = ±1. 
Once this embedding is done, our initial hamiltonian written in terms of continuous spins - here, hamiltonians (]?]) and 
(^) - can be rewritten as a generalized random-bond Ising model. In the Monte Carlo simulation the spin variables of 
this new problem are updated using an efficient non-local algorithm for Ising variables (e.g. the standard Swendsen- 
Wang algorithm). To flip the Ising variable corresponds to make a reflection of the vector Si in the hyperplane 
perpendicular to r. A necessary condition to get an efficient Wolff-type algorithm is that this transformation preserves 
the total energy of the system. In that case its application to a large set of spins costs only a surface energy and 
large-scale changes in the spin configuration are possible. Here, such a condition is verified only for the symmetric 
part (pi = ps) of the hamiltonians. As a consequence, we have chosen to perform our simulations with the reference 
symmetric hamiltonians. More precisely, for the 5*0(3) ® 0(2) model we consider: 

ij(°)=-^4 Pl ((^.5,) 2 -i) (35) 

<i,j> ^ ' 

and for the SU(2) ® 17(1) model we take 

# (0) = -J2 8 P^-S 3 (36) 

<i,3> 

Note that both hamiltonians reduce to the same hamiltonian in the spin-wave approximation (Hgyy — 
— 53<j SpiSSi.SSj where 6S represents the spin deviation from the reference vector). Note also that the zero- 
temperature ground-state energies are different for the two hamiltonians. This is not important since the various 
quantities computed in this work do not depend on this reference energy. 

Calculation of exact properties associated with the full hamiltonians are done by reweighting appropriately the 
Monte Carlo averages. Let us write Q the average of an arbitrary function Q(S) of the spin configuration S 
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Q = I dSi... I dS M Q(S)e~ H ^/ T /Z (37) 



where M is the total number of spins considered, H(S) is the exact hamiltonian, and Z the partition function. We 
re-express Q as follows 

Q = [<&!... [ dS M Q{S)e- v ^l T e- H{a) WI T l [dS x ...( dS M e- v ^ T e- H( ° } ^ T (38) 



where V is the difference between the exact and reference hamiltonian, V = H — H^ ' The Monte Carlo averages 
are performed over the set of spin configurations distributed according to the Boltzmann weight associated with the 
reference hamiltonian (denoted here as ((• • • )) ) 

Q = ((Q(S)e- v W' T )) /{(e- v W T )) . (39) 

By using such a procedure we are sure that at the symmetric point (pi = p^) our simulations are free from the 
critical slowing down problem (see, Ref.( fl32]| )). Away from the symmetric point, the situation is less clear. There is 
a subbtle interplay between the loss due to the undesirable fluctuations of the weight in averages (a growing source 
as the asymmetric parameter Ap = pi — pz is increased) and the gain obtained from the treatment of the large-scale 
collective-mode moves issued from the reference hamiltonian. In practice, we have found that realistic calculations can 
be done only for a small value of Ap with a maximum value of about 0.2. Although we have not made a systematic 
study of the effective (associated with H and not H 1 - ^) dynamical exponent, we are quite confident that, in the regime 
Ap < 0.2, the convergence of the estimators at the very low temperatures we have considered is very satisfactory. 
In order to compute the spin-stiffnesses we have used formulas (|l2|), ( |l8| ) and (^) presented in the previous section. 



V. MONTE CARLO RESULTS 



In this section we present the calculations of the spin-stiffnesses for the 50(3) (g> 0(2) and 5/7(2) <E> U(l) models. 
Before that, let us first give a qualitative visualization of the effect of topology in this problem. Figure (|l]) presents 
the specific heat as a function of the temperature for the two models and for three different lattice sizes: 4x4, 10x10 
and 20x20. Already for these relatively small systems the effect of topology is striking. In the 50(3) ® 0(2) case, 
the specific heat curves show a marked maximum whereas it is not the case for the topologically trivial 5/7(2) ® U(l) 
model. This maximum is usually interpreted as the signature of the presence of topological excitations. Note that 
the location o£ the maximum provides a rough estimate of the cross-over temperature at which these excitations 
are activated H . Here, and without making a detailed analysis based on much bigger sizes, we get approximately 
T v ~ 2.6. 
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FIG. 1. Specific heat as a function of the temperature for the SO(3) ® 0(2) and SU(2) <g> (7(1) models. Three different sizes 
have been considered: 4x4 (crosses), 10x10 (solid triangles), and 20x20 (solid squares). For the SU(2) Cg>t/(1) case the curves 
corresponding to the two largest sizes are undistinguishable. 



From our numerical results we propose to distinguish four different temperature regimes. 



A. The spin- wave regime 

The spin-wave regime corresponds to the very-low temperature regime. We have plotted in Figures (§) and (§) the 
spin-stiffnesses pi/T and p^jT as functions of lni at temperature T = 0.5 for the SO(3) <8) 0(2) and SU(2) <g> U(l) 
models, respectively. The parameters of the action are (pi,pi,P3) = (1,1,0.9). The behaviors of the spin-stiffnesses 
of both models are in full agreement with the two-loop RG predictions (solid line), Eq.(|3l|). As expected from our 
estimate of TV, these results show that at T = 0.5 the topological excitations are not yet activated and that the 
physics is controlled by spin-waves well described by the perturbative NLtr model. For this temperature and for the 
parameters pi and p% chosen, we find an almost linear behavior as a function of the logarithm of the size. In other 
words, the two-loop effects are almost negli gibl e. The numerical slopes are fully compatible with the theoretical slope 
of 1/4-7T as given by the one- loop equation (^3|). Note that, by measuring the two independent spin-stiffnesses of the 
HAFT model, Southern and Young H had already confirmed, for one value of the temperature, the RG predictions. 
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FIG. 2. 5*0(3) ® 0(2) spin-stiffnesses pi /T (lower curve) and pz /T (upper curve) as a function of ln L. Temperature T = 0.5 
and (pi,pi,ps) =(1,1,0.9). The solid line is the two-loop RG prediction as given by Eq. j3l|) . Open boundary conditions. 
Number of clusters used ranges from 5 10 6 to 8 10 6 . 
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FIG. 3. 5*17(2) ® i7(l) spin-stiffnesses pi/T (lower curve) and ps/T (upper curve) as a function of ln L. Temperature T = 0.5 
and (f»i,pi, P3) =(1,1,0.9). The solid line is the two-loop RG prediction as given by Eq.(j3l|. Open boundary conditions. 
Number of clusters ranges from 2 10 6 to 4 10 6 . 
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This overall behavior persists up to temperatures of order T ~1 where we enter a new regime. Finally, note that 
the absolute values of the temperature-rescaled spin-stiffnesses for the two models at T = 0.5 are different. Between 
T = and T ~ 1 this difference is almost constant. At T = 0, this constant can be calculated analytically, on the 
lattice models, from the finite parts of the one-loop counterterms that renormalize the couplings. It is found to be 
equal to 5/16. 

B. The almost-spin-wave regime 

For T ~ 1, the spin-stiffnesses of the two models start to differ: whereas their variations as a function of InL arc 
correctly described by Eq.([3l|) (see Fig.([|)) the absolute values of pi/T and pz/T for the 5*0(3) ® 0(2) model get 
smaller and smaller compared to those of the SU(2) ® U(l) model as the temperature increases. We have decided to 
refer to this regime as the almost-spin-wave regime. 
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FIG. 4. Spin-stiffness as a function of InL at T = 1.5 for the two models in the symmetric case (pi — ps). Note that in this 
case all three spin-stiffnesses are equal. 



Let us now give a quantitative account of this phenomenon. Since the difference between pi/T and pz/T turns 
out to be irrelevant for this discussion, we restrict ourselves from now on to the fully symmetric 50(3) (8> 50(3) and 
SU(2) ®SU{2) models where all three spin-stiffnesses are equal. The important point is that in the almost-spin-wave 
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regime, the spin-stiffness still displays a linear behavior as the function of lni even for the topologically non-trivial 
model. It is thus natural to define a characteristic length £ e // as: 



(40) 



with l; e ff being a function of T only. We recall that in the topologically trivial SU(2) ® SU(2) model and at the 
one-loop approximation, £ e ff is the correlation length (see Eq.(p3|)). Thus, the simplest hypothesis is that the £ e tf 
of the SO(3) <g> SO(3) model is still, in this regime, proportional to the correlation length. One can thus expect the 
ratio: 



R MCr T) = teff[SU(2)®SU(2)} 



(41) 



f e// [SO(3)®SO(3)] 

to be a good indicator of the influence of the topology. We give in FigE, R MC (T) as a function of the temperature 
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FIG. 5. Ratio of correlation lengths as a function of the temperature as denned in the text. The solid line is just to guide 
the eyes. The separation at T ~ 1 between the pure spin-wave (SW) and the almost-spin-wave regimes is rather arbitrary. 



From this figure, we see that at low temperatures (in the spin- wave regime) R M (T) is almost independent on 
the temperature. Note that at T = it converges to a value different from one because of the constant shift of 
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5/16 between the two temperature-rescaled spin-stiffnesses as discussed above. Within the range of temperatures 
T G [1, 2.1] the ratio is found to increase extremely rapidly. We have found that the curve can be well fitted using a 
form: 



R MC (T) = Cexp 



(T - T C )P 



(42) 



The "best" values found are: C = 23.52, a = 1.655, (3 = 0.600, and T c = 2.532. This clearly indicates that the 
topologically non-trivial model becomes more disordered than the topologically trivial one at a temperature of order 
T c = 2.5, a value compatible with that obtained from the peak of the specific heat. Of course, the analytical form 
chosen in Eq.(|4|) must be taken with lot of caution. Many different analytical forms could have been used and give 
similar results. 



C. The vortex regime 



At temperatures higher than typically T = 2.1, we enter in a new regime that we propose to call the vortex regime 
in which the topological excitations play a major role. Figures ||, 0, and ^ present, for the two models, the spin- 
stiffness as a function of the size at the temperatures T = 2.2,2.4 and 2.6, respectively. For the SU(2) ® SU(2) 
model (figure ||) the spin-stiffness is still correctly described by the two-loop NLc model predictions, Eq. (|3l|) . For 
the SO(3) Cg) 5*0(3) model it is still possible to define an effective slope at sufficiently small sizes but it now differs 
significantly from the RG predictions. This effective slope is found to increase as a function of the temperature. For 
example, at T = 2.2 the slope is about ~ —0.115 to be compared with —l/4ir ~ —0.080. At T = 2.4 the slope is 
~ —0.158, a value approximately two times larger than in the spin- wave regime. Moreover, the spin-stiffness displays 
some irregularities which could be associated with the presence of long-lived topological configurations that affect the 
dynamics. The most irregular curve has been obtained at temperature T = 2.6 (figure |^). 
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FIG. 6. 50(3) ® SO (3) spin-stiffness as a function of InL. T = 2.2. The solid line is the two-loop RG prediction. 
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FIG. 8. 50(3) <g> SO(3) spin-stiffness as a function of InL. T = 2.6 
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FIG. 9. SU{2) ® 5(7(2) spin-stiffnesses as a function of rnPfor^T = 2.2, T = 2.4 and T = 2.6. The solid lines are the 
two-loop RG predictions. 



D. The high-temperature regime 

At temperatures higher than typically T — 2.6 the spin-stiffnesses as a function of the size of the 5*0(3) ® 5*0(3) 
model recovers a smooth behavior. Fig. (|To|) presents such a behavior at temperature, T = 2.9. Fig.(pT|) presents the 
spin-stiffness for the SU(2) case at T = 7.7. The overall behavior of the spin-stiffnesses is quite different. In the 
50(3) ® 50(3) case p/T vanishes abruptly with a change of concavity whereas in the SU(2) (8 SU(2) case it goes 
slowly down to zero without any change of concavity. 
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FIG. 10. SO (3) ® 5*0(3) spin-stiffness as a function of InL. T — 2.9. High-temperature regime for the model. 
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FIG. 11. SU(2) ® SU(2) spin-stiffness as a function of InL. T = 7.7. High-temperature regime for the model. 



VI. DISCUSSION 



We now attempt to give a theoretical analysis of our results. It is notoriously difficult to tackle with the physics 
of Z2 topological defects since in this case there is nothing equivalent to the Villain transformation. In the usual 
ferromagnetic case, 0(N)/0(N — 1), Cardy and Hamber have proposed to describe, the effect of compacity of the 
sphere Sn-i by means of additional terms in the RG equation for the temperature EJ. These equations, valid in the 
vicinity oi N = 2 and at order T 3 - which correspond to two-loop in perturbation theory - read: 



dT{l) 
dl 



(N-2) 



T{lf 



2n 



(N 



T(l) 3 
W 



47r 3 y(0 2 + 



dy{lf / _ 2tt 



(43) 



dl 



T(l) 



y(i) 2 



These equations have been derived by assuming analyticity in y 2 and in N and by requiring that the two following 
limiting cases are recovered: i) the perturbative (3 function of the O(N) model for N > 3; ii) the Kosterlitz-Thouless 
equations for N = 2. In fact, Cardy and Hamber have shown that it is the only set of equations compatible with 
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these requirements H. For N = 2, one recovers the Kosterlitz-Thouless (K-T) equations Eil-Eil where y identifies 
with the fugacity of vortices. The case y = leads to the two-loop perturbative (3 function of the 0{N)/0(N — 1) 
NLcr model in two dimensions. For N ^ 2, y lacks of a clear interpretation but it has been conjectured that it encodes 
the effect of compacity. Clearly, the (N > 3, y ^ 0)-case we consider here could very well lie outside the domain 
of validity of these equations. It is therefore important to insist on the fact that they must only be considered as 
some phenomcnological RG equations for a model displaying topological defects, y playing the role of a fugacity by 
analogy with the K-T case. Our aim is to show that they are able to reproduce the gross features of the behavior of 
the spin-stiffness and correlation length found in our case. As in the K-T case, the physical "fugacity" y(l = 0) that 
appears as the initial condition in Eq.([l3]) is not independent on the temperature. However, in contrast with this last 
case its dependence on the physical temperature T(l = 0) is unknown. The simplest assumption we can think of is: 

y(l = 0) = e -^ T ('=°) (44) 

as in the K-T case. In Eq.(^), 7 is an adjustable parameter. 

The equations (|4|) are considered in the case N = 4 since 50(3) ® SO(3)/SO(3) and SO(4)/SO(3) have the same 
spin-wave content and differ by their topological properties. These latter properties are expected to be taken into 
account via the y-terms in Eq.(j43|). Moreover, we write them in terms of the spin-stiffness of the SO(3)®SO(3)/SO(3) 
model. These equations are obtained by the substitution 1/T(l) — > ip(l): 



' dp(l) _ 1 1 



dl 4tt 32tt 2 p(Z) 



mnhj(l) 2 p(l) 2 

(45) 
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By direct integration of Eqs.(|45|) up to the scale L of the lattice size - I =lni/a - one can obtain the dependence of 
p(l) on the temperature T(l = 0) and on L. To make contact with our Monte Carlo results, we are also interested in 
the correlation length £. This last quantity is defined, as usual, from the fact that when the RG scale e l becomes of 
the order of £, the spin-stiffness p(l) vanishes, see Eq. j4C|) : 

£/a ~ e l with p(l) ~ . (46) 

As in the Monte Carlo simulations, we compute from Eqs.(fl5|) the spin-stiffnesses and correlation lengths in both 
situations: without and with defects. This consists in setting respectively y(l) to zero or not in Eq.(|45|). Note finally 
that, for simplicity, we have chosen to take the same normalization at zero-temperature for the spin-stiffnesses with 
and without defects. 

Let us now show that we retrieve the essential features of the regimes previously identified except, obviously, for the 
high-temperature regime which is out of reach of the RG equations (fl5]) which, for the spin- wave part, are perturbative 
in the temperature. At very low temperatures, y(l) is very small and remains small along the RG flow. As a result we 
find almost no difference between the systems with and without vortices: this is the spin-wave regime that we recover 
trivially. 

As the temperature increases, the y-tenm plays a more and more important role. Adjusting the free parameter 7 
at a value 7 = 0.45 the temperature where the defects start to play a significant role is typically T(l = 0) ~ 0.15. We 
plot in Fig.fll^) the spin-stiffnesses as a function of the system size at a slightly higher temperature, T(l = 0) = 0.2. 
The upper curve corresponds to the system without defects and the lower one to that with defects. 
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FIG. 12. Spin-stiffness as a function of InL with and without defects. T(l = 0) = 0.2. 



As in our Monte Carlo simulations, we find that the absolute value of the spin-stiffness is decreased by the presence 
of defects. After an abrupt jump at very small sizes the spin-stiffness is found to be linear as a function In L. Up to 
an accuracy of a few percent, the slope is not affected by the defects. Its value, —0.080 ~ —1/Air, corresponds to the 
perturbative RG result. This behavior is similar to that predicted by the spin-wave analysis, except that the absolute 
value of the spin-stiffness is smaller. This corresponds to the almost-spin-wave regime previously identified. 

To obtain a completely consistent picture, it is also necessary to see whether the ratio of the correlation lengths, 
without and with defects, considered as a function of T(l = 0): 



R(T) 



€y(i)=o 



(47) 



behaves as R MC , Fig.(||). We give in Fig.fll3|), the ratio R obtained by direct integration of Eq. (j45|) . 
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FIG. 13. Ratio R(T) of correlation lengths as a function of T(l = 0). 



This figure clearly shows a behavior comparable to that observed in Fig.(|5|). 

At higher temperature, one enters in a regime where the spin-stiffness as a function of the size begins to display 
a different behavior. Figure ( |l4| ) presents the spin-stiffnesses at T — 0.5 (to facilitate the comparaison the dashed 
line represents the linear behavior associated with the model without defects) . It is possible to define a linear regime 
but now with a slightly greater slope than in the defect-free case. The difference between the two slopes is about 
10%. As the temperature is increased the spin-stiffness still displays a linear behavior but now with an increasing 
temperature-dependent slope. 
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FIG. 14. Spin-stiffness as a function of lnL with and without defects, T(l = 0) = 0.5. The dashed line represents the linear 
behavior associated with the model without defects. A slight increase of the slope is observed when the defects are present. 



Such a behavior has been numerically observed in the vortex regime (see, Fig.(||) and (0) where some "effective" 
linear behavior with a larger slope is observed). However, note that the irregular behavior found in Fig.(^) is not 
reproduced here. Of course, the range of temperature over which this linear behavior is observed, as well as the 
variation of the slope with T{1 = 0) depends rather strongly in our calculation on the relation between y(l = 0) and 
T(l = 0), Eq. (|44|) . With our choice of 7, the maximum variation found for the slope as a function of the temperature 
is about 15% which is somewhat below what is obtained in the simulations. This could certainly be corrected by 
another choice of y(T). 

At even higher temperature, Eqs.([45|) are no longer valid since they are based on a low-temperature expansion and 
the comparison with the numerical results does not make sense. 

The preceding analysis shows that the simple set of equations ([IB]) together with the relation (Q) seem to capture 
some of the important features of the presence of topological defects in the almost-spin-wave and vortex regimes. Of 
course, only a microscopical approach of the problem could allow to go beyond this semi-quantitative description. 
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